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A Cartesian adaptive level set method 
for two-phase flows 

By F. Ham and Y.-N. Young 


1. Motivation and objectives 

Simulations of flows involving free surfaces have become ubiquitous in the literature. 
The evolution of both simulation methods and computer speed have allowed the inves- 
tigation of many problems of increasing complexity and engineering relevance. For 2D 
simulations, it is reasonably common to use grid resolutions of 512 x 512 or larger. See 
for example the 2D planar breaking wave simulations of Chen et al. (1999), or the ax- 
isymmetric drop breakup simulations of Han and Tryggvason (1999a, 1999b). Many such 
free surface problems are fundamentally three dimensional and the absence of three- 
dimensional effects such as the pinch-off of stretched filaments by surface tension makes 
it difficult to make valid physical conclusions or use 2D results to develop models. Achiev- 
ing similar grid resolution in 3D, however, is outside the computational reach of most 
research centers. 

For many problems of interest, a fine grid is required in only a relatively small portion 
of the domain (e.g. around regions of high surface curvature) and much coarser grids are 
sufficient to capture the smoothly evolving velocity fields far from the surface. It is this 
type of problem that motivates the present Cartesian adaptive method, where the finest 
resolution in the neighbourhood of the free surface may be of order 512 3 (sometimes 
called the effective resolution), but the coarsening away from the free surface yields an 
actual number of control volumes that is less than using the fine grid everywhere by at 
least an order of magnitude. 

A number of adaptive Cartesian methods for freesurface calculations have been pro- 
posed in the literature. For example, Haj-Hariri et al. (1997) used Cartesian local grid 
refinement to investigate thermocapillary motion of 3D drops. Their method used an 
octree-based structure to store the grid heirarchy. More recently Sussman et al. (1999) 
proposed an adaptive levelset method based on the embedded grid approach, where the 
global computational domain is divided into uniformly refined rectangular patches. This 
approach has the advantage of allowing a traditional staggered or semi-staggered struc- 
tured discretization throughout most of the domain, with special interpolations required 
only at the interface between grid levels. The embedded grid approach may not, however, 
lead to the smallest Cartesian grids possible because of the restriction that refined regions 
be rectangular and isotropic. Additionally, embedded grid approaches can be difficult to 
parallelize and load balance. 

In the present contribution we develop a level set method based on local anisotropic 
Cartesian adaptation as described in Ham et al. (2002). Such an approach should allow 
for the smallest possible Cartesian grid capable of resolving a given flow. The remainder 
of the paper is organized as follows. In section 2 the level set formulation for free surface 
calculations is presented and its strengths and weaknesses relative to the other free surface 
methods reviewed. In section 3 the collocated numerical method is described. In section 
4 the method is validated by solving the 2D and 3D drop oscilation problem. In section 
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5 we present some results from more complex cases including the 3D drop breakup in an 
impulsively accelerated free stream, and the 3D immiscible Rayleigh-Taylor instability. 
Conclusions are given in section 6. 


2. Mathematical Formulation 

For recent reviews of the level set method and its many applications, the interested 
reader is referred to Osher and Fedkiw (2001) or Sethian (2001). Here we briefly present 
its formulation for incompressible 2-fluid problems. 

The incompressible, immiscible, two-fluid system is treated as a single fluid with strong 
variations in density and viscosity in the neighborhood of the interface. The continuity 
and momentum equations for such a flow can be written in conservative form as: 



dpui 

dt 


dpuj_ = 
dxj 
dpujUi 
+ dxj 


dp 

dxi 


dr Vy 

dxj 


+ pgt + an5{d)ni , 


(2.1) 

(2.2) 


where Ui is the fluid velocity, p the fluid density, p the pressure, the viscous stress 
tensor, gi the acceleration due to gravity, er the surface tension coefficient, k the local free 
surface curvature, 5 the Dirac delta function evaluated based on d the normal distance 
to the surface, and n, the unit normal at the free surface. 

It is well known that solving equations 2.1 and 2.2 directly in the presence of the 
discontinuities at the interface will lead to the excessive smearing of the interface over 
long time integration, and/or numerical “ringing” in the region of the discontinuities due 
to dispersive errors in the numerical discretization. In the present work, we capture the 
interface as the zero-level isosurface of a higher-order function </>, the level set function 
(effectively an additional transported scalar). On either side of the interface, the sign of 
the level set function can be used to determine which fluid is present. This allows the 
fluid properties to be written as a function of the level set only, specifically p = p(</>). 
This allows the continuity equation to be written: 
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= 0 


For the case of constant density except in the neighborhood 
which represents the fluid interface, the solution of the continuity 
decomposed into the solution of the following: 


(2.3) 

of the zero level set, 
equation can thus be 


dcj) 

dt 

dug 
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dcj) 

1 dxj 


= 0 




= 0 


(2.4) 

(2.5) 


Eq. 2.4 is the standard evolution equation for the level set function, and eq. 2.5 is the 
incompressible continuity equation. Note that the level set equation need only be solved 
near the interface, which we have defined as the zero level set. Away from the interface, 
it is common to make (f> equal to a signed distance function. 

In a similar way, the application of chain rule to the momentum equation allows the 
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Figure 1. Spatial location of variables for collocated discretization. 

decomposition of the momentum equation into the solution of the level set equation, eq. 
2.4, and the solution of the following momentum equation: 


Equations 2.4 - 2.6 comprise our governing system for the present work. Note that 
the level set formulation results in a system of governing equations that can no longer 
be written in conservative form. Thus, when the finite volume method is applied to this 
system, we cannot expect to achieve discrete conservation of mass and momentum (in 
the region of the interface). It is the hope of the level set formulation that these errors 
in conservation are mitigated by the more accurate capturing of the interface that is 
possible with the smoothly-varying level set function. 


3. Numerical Method 

The governing equations are solved on a Cartesian adaptive grid. The grid arrangement 
and adaptation are described more completely in Ham et al. (2002), and here we only 
discuss details specific to the simulation of multiphase flows using level set methods on 
such a grid. 

Figure 1 shows the spatial arrangement of variables on the Cartesian adaptive grid. 
All variables are stored at the control volume (CV) centers with the exception of a 
face-normal velocity, located at the face centers, and used to enforce the divergence- 
free constraint in each time step. The variables are staggered in time so that they are 
located most conveniently for the time advancement scheme. Denoting the time level by 
a superscript index, the velocities are located at time level t n and t n+1 , and pressure, 
density, viscosity, and the level set at time levels t n ~ s and t n+ 2 . The semi-descretization 
of the governing equations in each time step is then as follows. 

Step 1. Advance the level set: 





At 


(3.1) 




230 


Ham & Young 

Here the advantage of time-staggering is evident - i.e. we know the velocity w” from 
the previous time step, so this equation is decoupled from the other equations, and can 
be advanced on its own to get the level set value at the mid-point of the new time level, 
(j> n+ 2 . Here we have used an implicit Crank-Nicholson time advancement scheme, which 
requires an iterative solution method. In practice the hyperbolic system is not stiff, and 
can be quickly converged by a simple iterative scheme such as Gauss-Seidel. The spatial 
derivatives required in eq. 3.1 are approximated by a 5 tft- -order WENO scheme (Jiang & 
Peng 2000) . While effective at recovering the viscosity solution to the level set equation, 
the WENO discretization requires that the grid be locally refined in a uniform way in a 
band about the zero level set, and thus does not take full advantage of the flexibility in 
the adaptive grid. 

Step 2. Reinitialize the level set by solving to steady state: 

= sgn{4>) (1 - |V0|) . (3.2) 

The reinitialization step is performed every time step, and ensures that the level set is a 
signed distance function away from <p = 0, without (theoretically) changing the location of 
the zero-level set. Spatial derivatives required in eq. 3.2 are once again approximated using 
a 5 t/l -order WENO scheme (Jiang & Peng 2000). For this equation we use explicit third 
order TVD Runge-Kutta method for time integration. In practice, it is not neccessary 
to solve to steady state, and we have found that 5 iterations with a pseudo-time step of 
At = Aa:/4 is sufficient when performed every computational time step. 

Step 3. Calculate the new density and viscosity at the midpoint of the time step: 

p n+ * = piH (cj> n+i ) +Pi(\-H (3.3) 

p n+ * = Pi H (> +i ) + P2 (l - H (p n+i )) (3-4) 

where pi,pi are the constant properties of the fluid occupying the region (j> > 0, and 
p 2 , p 2 are the constant properties of the fluid occupying the region </> < 0, and H is a 
smoothed Heaviside step function (Sussman et al. 1994). 

Step 4. Project the momentum equations: 

The remaining steps are a variant of the collocated fractional step method described, 
for example, by Kim and Choi (2000). First, calculate a projected velocity field w* using 
the momentum equations: 


Uj - uf 

At 



(3.5) 


where Ri contains all other terms in the momentum equation, eq 2.6. We discretize both 
the convective and viscous terms implicitly using second-order symmetric discretizations, 
and the surface tension explicitly based on the known level set at the midpoint of the 
current time step, (f> n+ 2 . Some care must be taken in the discretization of the surface 
tension terms when treated as source terms in the momentum equations of a collocated 
scheme, as described below. 
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Step 5. Subtract the old pressure gradient: 


♦ Tl~\~ 1 
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1 dp 


' p n+ i Sxi 

Step 6. Interpolate the starred velocity field to the faces: 
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where () is a second-order interpolation operator that yields a face-normal component 
from two CV-centered vectors. 


Step 7. Solve the variable coefficient Poisson equation for the new pressure: 


J_ \ ' n *n+l A . 1 

At ^ 7 f n+ 1 

f f Pf 


n+i 


A 


/• 


(3.8) 


Step 8. Update the face velocities to the new divergence-free field: 


Un+l _ jjj 

At 


1 dp 
"+3 dn 

Pf 


n+ = 


(3.9) 


Step 9. Reconstruct the pressure gradient at the CV centers: 


_j_dp_ n+1 * _ 1 dp n+1 i x 

p n +h dxi dn 

r Ur 


(3.10) 


where () ' represents a reconstruction operator. In this case we use a face-area weighted 
least squares reconstruction. 


Step 10. Update the CV velocities: 


< +1 ~< 

At 


1 dp n+ i 

p n+ 5 dxi 


(3.11) 


By adding eqs. 3.5, 3.6, and 3.11, it is clear that a second-order time advancement of 
the momentum equation is recovered. 

A critical difference between the present formulation and the formulation of Kim and 
Choi is in the calculation of the starred face-normal velocities (eq. 3.7). Kim and Choi 
assume that: 


j R i "+ 1 / 2 / R n f + 1/2 

p n + 1/2 ~ n " +1 / 2 ’ 

r Ur 


(3.12) 
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This is an 0{ Ax 2 ) approximation, seemingly consistent with the overall accuracy of the 
method, and significantly simplifies the calculation of the Poisson equation source term. 
In the present investigation, however, it was found that when surface tension forces were 
introduced in the region of the zero level set, this approximation could lead to large non- 
physical oscillations in the CV-centered velocity field. To solve this problem, the surface 
tension forces must be calculated at the faces, and then averaged to the CV centers, i.e.: 



P Pf 


(3.13) 


With this calculation of the surface tension forces, we can no longer make the assump- 
tion of eq. 3.12, and the additional terms must be included in the calculation of UJ and 
thus in the Poisson equation source term. 


4. Validation 

As validation cases for the method, we solve the 2D column and 3D drop oscillation 
problems. 


4.1. 2D column oscillation 

From Lamb (1932), for a 2D column of liquid in the limit of zero viscosity perturbed 
according to: 


r = a (1 + ecos (nO)) (4.1) 

where a is the mean radius, e a small perturbation, and n the mode, the oscillation 
frequency will be: 


u> 2 = n (n 2 — l) (4.2) 

pa 6 

Figure 2(a) shows the initial condition for a drop oscillation calculation for mode 
n = 5. Here we have exaggerated the perturbation amplitude to e = 0.15 for illustrative 
purposes. In the actual computations, e = 0.02 was used for the initial condition. Figure 
2(b) compares the calculated results to the theory (period T = 2Tt/w), with generally 
good agreement. For these calculations, we used periodic boundary conditions in a square 
domain of size L = 4a, and set the density of the surrounding fluid to p 2 = O.OOlp to 
minimize its effect on the oscillations. The fluid viscosity was adjusted to keep the decay 
in amplitude of the oscillations to less that 5% per cycle. 

4.2. 3D drop oscillation 
For a 3D drop perturbed according to 


r = a (1 + eS n ) (4.3) 

where the surface harmonics of order n, S n , are defined: 



(4.4) 
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mode n 


(a) (b) 

Figure 2. a) Sample grid and initial zero-level set for 2D drop oscillation problem with n = 5. 
b) Comparison of computed and theoretical 2D drop oscillation periods for inodes n = 2 through 
5. 


Si = i vfi 3 ”* 2 "- 1 ) 

5 3 = - ( 5 cos^d — 3 cosO) 

54 = — (35 cos 4 9 — 30 cosd + 3) 

55 = -^7 \j^~ (63 cos 4 0 — 70 cos 2 9 + 15) 


(4.5) 

(4.6) 

(4.7) 

(4.8) 


the oscillation frequency is given by: 


lv 2 = n (n — 1) {n + 2) 


pa 


3 ' 


(4.9) 


Figure 3 (a) shows the drop surface for the mode n = 5. Figure 3 (b) compares the 
calculated results to the theory, also with good agreement. 


5. Results 

In the following subsections we present some results from simulations of more com- 
plex flows. For the purposes of this paper describing the numerical method, these cases 
are meant to simply illustrate the potential of the method. Consequently we do not 
present any analysis of the associated flow physics, which will be the subject of future 
investigations. 
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(a) 



(b) 


Figure 3. a) Initial zero-level set for 3D drop oscillation problem with n = 5. b) Comparison 
of computed and theoretical 3D drop oscillation periods for modes n — 2 through 5. 


5.1. Secondary breakup of spherical drops 

The secondary breakup of drops is an important problem in spray atomization, and the 
subject of the 2D axisymmetric numerical investigation of Han and Tryggvason (1999a, 
1999b). We used the present Cartesian adaptive method to calculate 3D drop breakup in 
an impulsively accelerated free stream. Figure 4 shows the calculated free surface for a 
case run with equivalent resolution of 128 x 128 x 128 with a domain size of 5 D x 5 D x 5 D. 

5.2. Rayleigh- Taylor Problem 

The Rayleigh-Taylor instability refers to the instability of the plane interface between 
two fluids of different density superimposed one over the other and subject to gravity. 
When the upper fluid has a density greater than the lower fluid, the interface can be 
unstable to small perturbations whose amplification is well described by linear theory 
with dependence on the density ratio, gravity, surface tension coefficient, and viscos- 
ity (Chandrasekhar 1961). As the amplification continues and the problem enters the 
nonlinear regime, the fluid mixing can become chaotic, characterized by bubbles of lighter 
fluid rising into the heavier fluid, and spikes of heavier fluid falling into the lighter fluid, 
with regions of high vorticity near the spike/bubble interface. The resulting mixing zone 
is known experimentally (Schneider et al. 1998) to broaden in a way that depends lin- 
early on the gravitational acceleration g and the Atwood number A = {p\ — P 2 ) / {p\ + P 2 ) , 
and quadratically on the time, i.e. 

hb,s = &b, s Agt 2 , (5.1) 

where hb or h s represents the penetration length for bubbles or spikes respectively, and 
a has been introduced as a constant of proportionality, sometimes called the acceleration 
constant. 

Recently, the miscible Rayleigh-Taylor problem has been investigated numerically by 
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(a) 0.5 (b) 1.5 (c) 2.5 (d) 3.5 (e) 4.5 



(f) 5.5 (g) 6.5 (h) 8.5 

Figure 4. Breakup of an initially spherical drop in an impulsively accelerated field. The 
number below each figure indicates the non-dimensional time tU/D. Free stream veloc- 
ity is from top-left-back to lower-right-front. Properties for this simulation are as follows: 
Re = pgDU/i-ig = 1000, We = p g U 2 D/a = 20, Oh — pi/C pm D = 0.007. 

Young et al. (2001), who calculate ab « 0.03. Experimental results for immiscible fluids 
yield higher values of ab ~ 0.05 — 0.07. A recent review of experimental and simulation 
results for the immiscible case is given by Glimm et al. (2001). 

Figure 5 presents some results from the application of the present Cartesian adaptive 
method to this problem. These computations were performed in a 2 x 2 x 2 box with 
the interface perturbed randomly with a maximum amplitude of 0.01. Properties of the 
simulation were selected to give a wavelength of maximum instability of A* = 0.35. These 
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Figure 5. Evolution of the interface for the 3D Rayleigh- Taylor problem with multi-mode 
perturbation: top panels - plan view looking down on interface; bottom panels - elevation view 
from side. 

preliminary simulations were made with the relatively coarse equivalent resolution in the 
region of the zero-level set of 64 x 64 x 64, or about 10 CVs per wavelength for the most 
unstable mode. Figure 5 illustrates the evolution of the interface at 3 different times for 
the Atwood number A = 0.98. 


6. Conclusions 

A simulation tool that integrates Cartesian adaptive grids with the level set method 
has been developed. The method as described is particularly suited to problems where 
the smallest scales are associated with the free surface motions, and significant savings 
can thus be realized by coarsening the grid away from the surface. The method has been 
validated and its potential demonstrated by solving several test problems, including 2D 
and 3D drop oscillations, drop breakup in an impulsively accelerated free stream, and 
the immiscible 3D Rayleigh-Taylor problem. 

Our ongoing work is focused on changing the level set solver to the particle level set 
method of Enright et al. (2002). The improved conservation properties of this method 
may allow us to move away from the higher-order WENO discretizations presently used 
in the level set advancement and reinitialization equations. With lower-order methods, 
adaptation can be added in the region of the zero-level set, resulting in a further signifi- 
cant reduction in computational cost. 
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